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Abstract 



This paper proposes a new method, in the frequency domain, to define absorbing 
boundary conditions for general two-dimensional problems. The main feature of the 
method is that it can obtain boundary conditions from the discretized equations 
without much knowledge of the analytical behavior of the solutions and is thus 
very general. It is based on the computation of waves in periodic structures and 
needs the dynamic stiffness matrix of only one period in the medium which can be 
obtained by standard finite element software. Boundary conditions at various orders 
of accuracy can be obtained in a simple way. This is then applied to study some 
examples for which analytical or numerical results are available. Good agreements 
between the present results and analytical solutions allow to check the efficiency 
and the accuracy of the proposed method. 

Key words: Absorbing boundary conditions, waveguide, finite element, periodic 
medium. 
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1 Introduction 



Wave problems in unbounded media can occur in many applications in me- 
chanics and engineering such as in acoustics, solid mechanics, electromag- 
netics, etc. It is well known that analytical solutions for such problems are 
available only for some special cases. On the contrary, numerical methods can 
be applied to many complex problems. Physically, for problems in infinite do- 
mains, the energy is produced by sources in the region to be analyzed and 
must escape to infinity. For methods solving the problem on a bounded do- 
main like the finite element method, it introduces the difficulty of an artificial 
boundary to get a bounded domain. This boundary must be such that the 
energy crosses it without reflection and special conditions must be specified at 
the artificial boundary to reproduce this phenomena. Generally, these can be 
classified into local or global boundary conditions. With a global condition all 
degrees of freedom (dofs) on the boundary are coupled while a local condition 
connects only neighboring dofs. 

The first global method which has been used for solving such problems was the 
boundary element method. This method is well adapted for infinite domains 
and is described in numerous classical textbooks like [1,2.3,4,5J. It consists 
in solving an equation on the boundary of the domain only and the radiation 
conditions are taken into account analytically. It also reduces the dimension of 
the problem to a surface in 3D and to a curve in 2D decreasing thus the size of 
the linear problem to solve. However, the final problem involves full matrices 
which are also generally non symmetrical. It is also mainly limited to lin- 
ear problems and to homogeneous domains or otherwise one has to introduce 
special and complex techniques to deal with non linear or non homogeneous 
situations. There are also singularities in the integrals which need special at- 
tention for the numerical integrations. So this method is interesting and has 
been extensively used but it can lead to heavy computations when the number 
of degrees of freedom increases. More information on such techniques can be 
found in the historical and review papers [6fT] . 




In the other approaches, the computational domain is truncated at some dis- 
tance and boundary conditions are imposed at this artificial boundary. These 
conditions at finite distance must simulate as closely as possible the exact 
radiation condition at infinity. An approach leading to a global boundary con- 
dition is the Dirichlet to Neumann (DtN) mapping proposed by [8|9] and in an 
earlier version by [TUlfTT] . It consists in dividing the domain into a finite part 
containing the sources and an infinite domain of simple shape. The solution 
in the infinite domain is solved analytically, for example by series expansions, 
and an exact impedance relation is obtained on the boundary between the 
finite and infinite domains. This relation links the variable and its normal 
derivative on the whole boundary. The DtN mapping is thus non local and 
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every node on the boundary is connected to all other nodes. This gives a full 
matrix for the nodes of the boundary which partially destroys the sparse ma- 
trix of the FEM and increases substantially the computing resources needed 
to get the solution. The solution has to be found in the exterior domain by 
analytical or numerical methods. When the analytical solution can be found, 
it is generally under the form of a series expansion. The number of terms in 
the expansion must be sufficient for an accurate solution which can lead to 
heavy computations. Developments of the method can be found in [T2][T3] . An 
application to the case of wave scattering in plates is also found in [T3] . 

The other methods are local and the condition at a node involves only neigh- 
boring nodes which make them less demanding in computing resources and 
much easier to implement in a finite element code but also less accurate. A 
first possibility of such approaches is the use of infinite elements as proposed 
by [Toll lb'lll7|l 1811 19j . It consists in developing special elements with a behavior 
at infinity reflecting that of analytical solutions obtained for the same type of 
problems. For wave problems, it involves complex-valued basis functions with 
outwarding propagation wave-like behavior in the radial direction. The ele- 
ments were further developed by [20] to considered other coordinate systems 
such as expansions in prolate coordinates. This method is interesting but the 
inclusion of infinite elements requires the development of special elements and 
these elements can depend on decay parameters which have to be accurately 
chosen. A review of these methods has been proposed by [2~T] . 

In the perfectly matched layer proposed by [22.23J, originally for electromag- 
netic waves, an exterior layer of finite thickness is introduced around the 
bounded domain. The absorption in this domain is increasing as we move 
towards the exterior such that outgoing waves are absorbed before reaching 
the exterior truncation boundary. The number of elements in the layer, its 
thickness, the variation of the absorption properties have to be carefully cho- 
sen to optimize the efficiency of the method. This efficiency is better for a layer 
with a large thickness but this can lead to a significant increase in the number 
of elements in the finite element model. Various developments of the method 
can be found in [24j and its optimization in [25]. Otherwise, various classes of 
absorbing boundary conditions were also developed by [20] • They consist in 
the numerical approximation of differential operators on the boundary. For in- 
stance, examples of the application of Bayliss-Turkel conditions are presented 
by [2?]. However, the more accurate boundary conditions involve high order 
derivatives on the boundary which are difficult to implement in the finite ele- 
ment method [2B]- Finally, it was proved by [25.29J that, in fact, the perfectly 
matched layer and the absorbing boundary conditions were closely connected. 
The Helmholtz equation was also solved with these two boundary conditions 
by [30] and these conditions were compared and optimized to minimize the 
reflection. Other boundary conditions involving only second order derivatives 
have also been proposed. They introduce auxiliary variables and systems of 
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equations on the boundary which lead to high order boundary conditions, see 
[31] for a review of such non-reflecting boundary (NRBC) methods. They were 
mainly developed for acoustic problems but in [32J a local boundary condition 
for elastic waves has been proposed. In [33] an impedance boundary condi- 
tion in new coordinates was developed for the convected Helmholtz equation. 
For fluid dynamic problems, [34] developed Lagrange multipliers for imposing 
various absorbing boundary conditions for cases where the type and the num- 
ber of boundary conditions can change, for instance as the flow changes from 
subsonic to supersonic regimes and its direction varies with time. A general 
review of the methods described in the precedent paragraphs for various dy- 
namic, acoustic and wave propagation problems can be found in [3511361137(38] . 
Comparisons are also made between the different methods. 

In the present study, another local method is proposed. This method works on 
discrete systems directly, in contrast with many existing absorbing boundary 
conditions which are written on the continuous differential equations and dis- 
cretized after. The principle of the method is to compute wave propagations 
in groups of elements near the boundary from the dynamic stiffness matrix 
of these elements. Then, a boundary condition is obtained for cancelling the 
reflected waves. This condition is finally written as an impedance boundary 
condition relating the force and displacement degrees of freedom on the bound- 
ary. The approach is based on the waveguide theory proposed by [39|40|41j and 
is used to determine absorbing boundary conditions at the truncation bound- 
ary of 2D periodic media. Only information related to one period, obtained 
from any standard FE software (the discrete stiffness and mass matrices and 
the nodal coordinates) are required to formulate the method. The advantage 
is that it can be applied to media with various complex behaviors. 

This paper is outlined as follows. In section 2, the methodology for determin- 
ing absorbing boundary conditions for periodic media is described. Then, a 
discussion for the application of the method to general media is proposed. In 
section 3, a simple application is proposed to show the results of the method 
in a case where detailed computations can be done. In section 4, two examples 
of finite element computations in acoustics and elastodynamics are presented. 
They allow to check the efficiency and accuracy of the proposed method for 
more complex cases. Finally, the paper is closed with some conclusions. 



2 Absorbing boundary conditions 

We suppose that we want to solve a mechanical problem on an infinite do- 
main exterior to the bounded domain Qi nt (see figured]). The infinite domain 
is approximated by the finite domain Q which is exterior to VLi nt and is limited 
by the exterior boundary T ext . We are looking for a solution with radiating 
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condition at infinity which means that the solution should be outgoing near 
the boundary T ext . Near this exterior boundary the solution can be seen as 
composed of incident waves denoted A + and reflected waves A_. For a per- 
fectly absorbing boundary, one should have = 0. In fact this condition is 
very difficult to implement in the numerical solutions of such problems. In- 
deed, only the global solution is easily computed but the decomposition into 
incident and reflected waves is difficult to obtain. The problem is thus to find 
an appropriate boundary condition to impose on the exterior boundary to 
finally get A_ « on the solution. To be easily included in a finite element 
model the searched boundary condition should be local and the condition at 
a node of the boundary should involve only neighboring nodes. 

The approach proposed in this paper consists in studying this problem by 
first considering the case of periodic media. For this case, positive and nega- 
tive waves and their amplitudes A + and A_ can be computed by the method 
presented below. Then an exact boundary condition can be formulated for a 
half-plane boundary. It is further shown how this condition can be approx- 
imated by a local condition on the boundary. As homogeneous media are 
special cases of periodic media, the method presented here applies also to ho- 
mogeneous media. Before considering the general simple example for 
the Klein Gordon equation will be presented. 



2.1 A simple example 



Consider first the stationary Klein Gordon equation given by 

3"4r + (k 2 - m 2 )u = 
ax 2 



where u is the solution and k, m are real parameters. This equation is dis- 
cretized with linear two nodes elements such that 

u(0=«lWl(0+U2#2(0 (2) 



where Aq(£) = (1 - f)/2, iV 2 (0 = (1 + O/ 2 , «i and u 2 are the values of the 
function at the two nodes of the element. The discretization of the first and 
second parts of relation ([1]) leads, for an element of length I, to the matrices 



1 -1 
-1 1 



m 



2 1 
1 2 



(3) 
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and the dynamic stiffness matrix of one element is given by 



1 -1 
-1 1 



o 



2 1 
1 2 



Waves of propagating constant are such that 
u 2 = e^u\ 

leading in an element to 

e 2 »di2 + (d u + d 22 )e^ + d 21 = 

where dij are the components of the matrix d. Taking into account the 
metries in the matrix d, this yields 

e 2» + 2 rfii e , + 1 = Q 
di2 

whose solutions are 



with 



dn l-/ 2 (fc 2 -m 2 )/3 
d^ 2 ~ ~l + l 2 (k 2 -m 2 )/6 

For l 2 (k 2 — m 2 ) « 1, one gets, at first order, 
e M± 1 ± i/V/c 2 - m 2 



meaning 



/i ± ps ±iWk 2 — m 2 
There are thus two waves in an element, such that 







1 


and 










d n + e^+rf 12 




./-_ 





tin + e^-rf 
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and the general solution is given by 



u 




u+ 








= a + 




+ a_ 




j_ 




J+. 







(14) 



The condition for only outgoing waves is thus a_ = on the right boundary 
and a + = on the left boundary leading respectively to the conditions 



f/u = f + /u + on the right 

jju — j-ju- on the left (15) 
The condition in the first case is 



f/u = d 11 + e> l+ d 



12 



d\2 V ^12 



2 

12 



'-(fc 2 -m 2 ) + — ((P -m 2 )/) 2 

_1_ _ 



= iVA; 2 - m 2 y 1 - ^(A; 2 - m 2 )/ 2 (16) 

miVk 2 — m 2 (17) 
In the second case, one gets 

//« = d u + e^-d 12 « -iVA; 2 - m 2 (18) 



We recognize approximations of the classical absorbing boundary conditions 
which have been obtained here directly from the discretized equations. Com- 
pared to the classical boundary condition on the right (and exact in this case) 
f/u — i\Jk 2 — m 2 , the relative error is ^1 — -^{k 2 — m 2 )l 2 which depends 
mainly on the size of the element relatively to the wavelength. The present 
boundary condition has been obtained entirely from the discrete matrices with- 
out any knowledge of the analytical solution of the problem. 

To estimate the reflection coefficient created by such a boundary condition 
conside r an i ncident wave Ae l ^ k2 ~ rn2x on the boundary. A reflected wave 
RAe~ l ^ k2 ~ m2x i s created. The total solution and its associated force are given 
by 
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u( y x)=Ae Wk2 - m ' 2x + RAe 

Writing the boundary condition (|T6|) . for instance by taking the boundary at 
x = 0, yields 



44 = - — ~ = \h~—{k 2 - m 2 )l 2 (20) 
u(0) 1 + R V 12 V ; K J 

So the reflection coefficient is finally given by 



r l-^l-j- 2 (k 2 -m 2 )P 



f + ^l-^(P-m 2 )/ 2 
l(fc 2 -m 2 )/ 2 (21) 



This coefficient is low and of second order when y/k 2 — m 2 l < < 1 . 
2.2 General impedance boundary condition 



In this section we present the general outline of the method before starting a 
more rigorous developement in the following section. So, to extend the prece- 
dent example to more general cases, consider a vector function u(x, y) and a 
force vector f(x, y) acting on a line parallel to the y axis as in figure El They 
can be decomposed by a Fourier transform as 



+oo +oo 

u(x,y) = J J u(k x ,k y )e^ x+ ^dk x dk y 



— OO — oo 

+0O +0O 



f(x,y) = 7? ^= I I i(k x ,k y )e^ x+k y^dk x dk y (22) 



(27T) 



-oo — oo 



Let us suppose that u is solution of a linear operator 



n=N o„ u 

L(u)=E £ a Qia2 ^— = (23) 

n=0 ai+a 2 =n ud/ u i> 



In the Fourier domain, this relation yields 

(n=N \ 
J2 E a *i*2 (ik x ) ai (^)° 2 ufo, fc,) = (24) 
n=0 ai+a2=n / 
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For a given value of k y , the precedent relation has non zero solutions for k x 
such that the determinant 



n=N 

E E 

n=0 ai+«2 = 



(ik x ) 1 {iky) 



a 2 



(25) 



Let us denote by kf the n + positive solutions such that Re(kj~) < or 
Re(kj~) = and the energy flux is directed towards positive values of x. We 
denote by kj the other solutions. We have the decomposition 

j=n+ j=n- 

u(x, k y ) = X «:u:<* * + £ ajuje*!* (26) 

j=l 3=1 



In the same way, for the force components 



j=n + j=n_ 

f (x, fcj,) = £ a+f+e^ ^ + £ a-f-e^ 7 * (27) 
i=i i=i 



where fj- and f~ are the force components respectively associated to and 
uj. If the boundary is such that only positive waves exists at proximity, one 
has 



3=n+ 

u(0, k y ) = £ a t u t = U + a+ 

3=1 
3=™+ 

f (0, ky) = J2 otf+ = F + a + (28) 
j'=i 

where U + and F + are the matrices whose columns are respectively u+ and 
fj~. Eliminating the a + coefficients, one gets 



f(0,k v ) = F + (U + )- 1 u(0,k v ) 

f(0,y) = (G*u)(0,y) (29) 

with 

G(ky)=F + (U + )- 1 (ky) (30) 



and * means the convolution. 



In the following this boundary condition will be computed directly from the 
discrete equations for general linear media. 
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2.3 Solution in a periodic medium 



Consider an infinite two dimensional periodic medium, as shown in figure [3j 
The elementary period is limited by the domain (a^a^) G [0,6i] x [0,62]- A 
function U(xi,x 2 ) defined on the two-dimensional periodic medium can be 
decomposed as an integral of pseudo periodic functions 

7T 

U(x 1 ,x 2 ) = J e ikX2 U(x h k,x 2 )dk (31) 

_7T_ 

i>2 



where U(xi, k, x 2 ) is a periodic function in x 2 with period b 2 . From the Fourier 
transform U(x\, k) of U(x\, x 2 ), one has 



U(xi,k,x 2 ) = — J2 U(xi,k + 2ir—)e 



27C m 2 =-oo ' b 2 



b 2 



rOO 



+OO 



1 \ l2K-r^-.r'' I -I I A" -|- 2 TT ^— I.I- , r , , 1 

— 2^ e ^ J a ^ U{xi,x)dx 



m,2=— oo 



+00 



-00 



— / V 5(x 2 - x - m 2 b 2 )e~ tkx U(xi,x)dx 

27F 7 m 2 =-oo 
+00 

£ e- ifc ( a:2+m262 V(a; 1 ,x 2 + m 2 & 2 ) (32) 



6 +v 



2tt 



m2=— 00 



This gives the relation inverse of fl3~T|) . From relation fl3Tl) . one sees that the 
behavior in x 2 of the general solution can be obtained from functions as 
e tkX2 U(x 1 , k,x 2 ) with U(xi,k,x 2 ) periodic in x 2 . Along direction 1, we use 
a decomposition in Bloch waves as it is usual in periodic media. Finally, the 
general solution can be obtained from functions u(x\, k, x 2 ) = e tkX2 U (xi, k, x 2 ) 
such that: 



u(x u k, x 2 + m 2 b 2 ) = e ifcm2b2 M(xi, k, x 2 ) (33) 
u(x 1 + mibi, k, x 2 ) = e imifl u(x 1 , k, x 2 ) (34) 



where m\ and m 2 are integers, k G M. H 



71 7T 

b 2 b 2 



and \i G C. 
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The discrete dynamic equation of a cell (an elementary period) obtained from 
a FE model at a frequency oj and for the time dependence e~ luJt is given by: 

(K - iuC - cu 2 M)q = f (35) 



where K, M and C are the stiffness, mass and damping matrices, respectively, 
f is the loading vector and q the vector of the degrees of freedom (dofs). 
Introducing the dynamic stiffness matrix D = K — iuC — co> 2 M, decomposing 
the dofs into boundary (B) and interior (I) dofs, and assuming that there are 
no external forces on the interior nodes, result in the following equation: 





















qz 








(36) 



The interior dofs can be eliminated using the second row of equation (|36|) . 
which results in 

q/ = -D7/D /B q B (37) 

The first row of equation (1361) becomes 

= (pBB - D B/ D7/D JB ) q B (38) 



which can be written as 

f = Dq (39) 



It should be noted that only boundary dofs are considered in the following. 

The periodic cell is assumed to be meshed with an equal number of nodes on 
their opposite sides. The boundary dofs are decomposed into left (L), right 
(R), bottom (B), top (T) dofs and associated corners (LB), (RB), (LT) and 
(RT) as shown in figure HI The longitudinal dofs vector is defined as q; = 
* [*qjr 'q^j t qLB t( iRB 1< \rt *qr/r]- Thus, equation (139^) is rewritten as 



D H 


B lB 


Bit 




q« 






Bbi 


Vbb 


Dbt 




qs 






D„ 


T>TB 


D^j 1 




q T 




f T 
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Using the pseudo periodic condition (1531) and the effort equilibrium at the 
bottom side of the cell, relations between the transverse dofs are given by 



q T = e ikb2 q B 



i B + e~ ikb H T = 



(41) 



Multiplying the third row of equation (1401 with e lkb2 , taking the sum of the 
second and third rows of equation (l4"01 . using conditions (jUJ), lead to 

(Dbj + e- lkb2 B Tl ) qi + (p BB + D TT + e~ ikb2 B TB + e ikb2 B BT ) q B = 0(42) 



so 



q B = - (D BB + D TT + e~ tkb2 T> TB + e lkb2 B BT ) 1 (D ra + e~ lkb2 Ti T i) q,(43) 



Using (|4T!) and (143]) . the first row of equation (1401 becomes 



D„ — ( D/ B + e lkb2 T> lT ) (D BB + D TT + e~ tkb2 B TB + e lkb2 B BT 



x 



(n Bl + e- lkb2 B Tl 



q/ 



(44) 



which can be written as 



(45) 



Using the pseudo periodic conditions (1341 also lead to the following relations 
between longitudinal dofs 



p i(fj,+kb2) 



(46) 



From the pseudo periodic conditions (T4"61 . it can be seen that all components 
of the vector qi depend on the set of dofs defined by q r = * [*qx *qL_e]- This 
can be expressed as 



q l= Wo + e^Wx q 



(47) 
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where the matrices Wo and Wi depend on the wavenumber k and are 
by 



I 


o 




o 


o 


o 


o 




I 


o 


o 


I 


Wi = 


o 


o 


o 


o 




o 


I 


o 


o 




o 




o 


e ikb 2 j 




o 


o 



The equilibrium conditions between adjacent cells are given by 

eH L + f R = 
e^LB + f RB + e^- kb ^f LT + e- ikb ^ RT = 

that can be written as 

(e ifi W* + Wt) fi = 

where (.)* denotes the operator of complex conjugate and transpose. 
Combining (jl7]) and (|5D|1 ; lead to 

(e^W* + W*) D, (W + e^Wi) q r = 

that can be written as 

(a + e^(Ai + A 2 ) + e 2 ^A 3 ) q r = 

where 

A = W^W 
Aj = W*D,W 
A 2 = W* 1 B l W 1 
A 3 = W*D ; W! 
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The eigenvalue e 1 ^ and the eigenvector q r are thus solutions of a quadratic 
eigenvalue problem. It is convenient to transform the problem ( |52|) into another 
linear eigenvalue problem as 





A 3 


O 




q r 






O 


A 3 _ 









O A 3 

-A -(Ai + A 2 ) 





q r 




q r 



(54) 



with q r = e iM q r . 

From equations (jHJ) and fj45l) . one can notice that 
'D,(fc) = D,(-*) 



(55) 



Moreover, from (148]) . we have 



for j 



0,1 



(56) 



and from (J53|) 



*A (/e) = A 3 (— fc) 
*A 1 (A;) = A 1 (-A;) 
*A 2 (A;) = A 2 (-A;) 
*A 3 (A;) = Ao(-A;) 



(57) 



It can be easily shown by taking the determinant of the matrix in relation (1521) 
that if e 1 ^ is an eigenvalue for the wavenumber k, e~* Mj is also an eigenvalue 
for the wavenumber —k. These represent a pair of positive and negative-going 
waves, respectively. The In eigensolutions of equation fl54l can be split into 
two sets of n + and n~ eigensolutions with 2n = n + +n~ , which are denoted by 
(e^ + ,q+) and (, 

In the case 



e 1 ^ 



< 1. 



l; , respectively, with the first set such that 

1, the first set of positive-going waves must contain waves 

propagating in the positive direction such that Re jicuq^fj j > where fj is 
the reduced set of boundary force dofs of left cells on right cells and is given 
by 



f 

3 



fiB + e- ifcfe2 f LT 



W*fJ 



WJD, Wo + e^Wx q^ 



(5f 
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In the second set of negative-going waves, the eigenvalues e 1 ^ are associated 
with waves such that Re jicuq^fj j < 0. 

With the eigenvector q, and the force component of relation ( 1581) . we introduce 
the state vector 



qj(fc) 










(Ai(fc) + e i ^^A 3 {k))q J {k) 



(59) 



In this relation cy(k) is the eigenvector associated to e w ®. One can also 
introduce 



y 3 {-k)= t p J (-k)(A 2 (k) + e i ^A 3 (k)) %{-k) 



(60) 



In this relation Pj(— k) is the eigenvector associated to e~ lflj ^ since we have 
seen that e~ % ^ 3 ^ is also an eigenvalue of (152]) for the wavenumber —k. From 
relation fl52|) written for the eigenvector cij(k), multiplying this relation by 
e -v*j(*0 anc i then on the left by t p i (—k), one gets 

* Pi (-fc) (e-^«A (A;) + (A^fc) + A 2 (fc)) + e^A^k)) qj(k) = (61) 

In the same way, writing relation (1521) for the eigenvector Pi(— k), taking the 
transpose of the relation, using relations fl57|) and multiplying on the right by 
qj(fc), leads, after a global multiplication by e % ^ k \ to the following relation 

'Pii-k) (e^A 3 {k) + (A 1 (k) + A 2 (k)) + e-^A (k)) ^(k) = (62) 



The difference between the two precedent relations yields 

( e <«(*) _ e w(*))tp.(_jfe) (a 3 (A;) - e -^ (fc) e-^ (fc) A (£;)) ^-(jfe) = (63) 



In the case e* Mi(7c ) 7^ e % ^^ k \ we get 

*Pi(-Jfe) (e^ (fc) A 3 (A;) - e^ (fc) A (£;)) c^(k) = 



(64) 



Now it is possible to compute the product Vj(— k).Xj(k) by 



y i {-k).x j (k) = t p i (-k){A 2 (k) + e i ^A 3 {k))^ 
+ t Pi(-A0(A 1 (fc) + e^A 3 (k))q 3 (k) 



16 



= * Pi (-fc)(A 2 (fc) + e^Az(k))^{k) 

-* Pi (-A;)(A 2 (A;) + e~ 1 ^ A {k))^{k) 
= didij (65) 

The result of relation flB^l) has been used in the case e^^' ^ e Vj '^ and c?j is 
a factor depending on the eigenvector %. This gives orthogonality relations on 
the statevectors associated to the eigenvalues. 



2.4 Absorbing boundary conditions 



Figure [2] presents the periodic medium near the exterior boundary. In this 
domain the solution is described by relation (131j) . yielding, respectively for 
the displacement and force components, 



&2 

q(xi,x 2 )= J q(x 1 ,k,x 2 )e lkX2 dk 

— JL. 

7T 

f(x u x 2 )= J f(x 1 ,k,x 2 )e ikX2 dk (66) 



b 2 



with the force components given by relation fl58|) . Introducing the state vector 
x = *(*q, *f) and decomposing this solution into the different waves, we get 




The last relation is the approximation obtained by the finite element com- 
putation of wave solutions presented before. The condition of outgoing waves 
means that there is no incoming wave, so the amplitudes cij{x\, k) associated 
with incoming waves must equal zero. This condition is obtained by 

j=2n 

yT(- k )- a A x ^ k )*j( k ) =° for l<l <n~ (68) 
3=1 
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In this relation yf(—k) are the vectors associated to the negative going waves, 
given by relation (1601) . Using relation (1651) . one gets a~(xi, k) = for 1 < j < 
rT for the amplitudes of the negative going waves. Introducing the matrix Y 
with lines given by y ; ~ leads to 

Y(-k).x(x 1 ,k,x 2 ) = (69) 



Decomposing now x into its displacement and force components, doing the 
same thing for Y(— k) = [Q(— k) F(— k)] leads to 

Q(-A;).q(x 1 , k, x 2 ) + F(-k).f(xi, k, x 2 ) = (70) 



The relation on the boundary is 

f ( Xl , k, x 2 ) = -F- 1 (-A;)Q(-A ; )q(x 1 , k, x 2 ) (71) 

and then from relation (166]) 

K 

f(xi,x 2 ) = - J F- 1 (-A;)Q(-fc)q(a; 1 ,fc,a; 2 )e^ 2 ^ (72) 

_7T_ 

From the inverse relation ( l32l) . one also has 

h +0 ° 

${x 1 ,k,x 2 ) = 7 r £ e- ik ^ +m ^q(x 1 ,x 2 + m 2 b 2 ) (73) 



which leads to 



m2=— oo 



6 9 "9 +°° 



f(xi,x 2 ) = ~ f F-\-k)Q{-k) Yl e- lk ^ +m ^q{x 1 ,x 2 + m 2 b 2 )dk{74) 

^ 7f _j L m 2 =-oo 

Introducing the function 

G(x 2 )=-7^ / F- 1 (-A;)Q(-A;)e-^ 2 ^ (75) 

Z7T J 



The final relation is 

+00 

f(x 1 ,x 2 ) = G{x 2 + m 2 b 2 )q(x 1 , x 2 + m 2 b 2 ) (76) 
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This is the impedance relation on the boundary obtained with the assumption 
that there is no negative going wave. This relation is the absorbing boundary 
condition we were looking for. It can be computed from the wave vectors and 
the force components associated with them. Relation (I75|) involves an infinite 
number of terms on the boundary. This relation can also be written as 



+00 



f(x l7 x 2 )={ Y G(x 2 +m 2 b 2 ) j q(:ci,x 2 ) 

\m2=-oo / 
+00 

+ Y G ( x 2 + m 2 b 2 ) (q(xi,x 2 + m 2 b 2 ) - q(x 1 ,x 2 )) 



7712 = — OO 
f +OO \ 

Y G(x 2 + m 2 b 2 ) J q(xi,x 2 ) 

Vm2=— 00 / 
J / +00 



+ 2b~\ E m 2b 2 G(x 2 + m 2 b 2 )] (q(x 1 ,x 2 + b 2 ) - q(x 1 ,x 2 -b 2 )) 



\ 7712 = — OO 



-OO 



+ Y G(x 2 + m 2 b 2 )[q(x 1 ,x 2 + m 2 b 2 ) - q(xi,x 2 ) 

7712= — OO 

-^m 2 (q(xi, x 2 + 62) - q(xi, x 2 - b 2 ))} (77) 

If q is slowly varying the last term should be small and for practical purposes 
we will use the approximate relations at various orders given by 



Q 

f{x 1 ,x 2 ) « G q(a;i, x 2 ) + — -(q(xi, x 2 + 62) - q(xi, x 2 - b 2 )) 

Z0 2 



G 



+7T72" (q(zi, ^2 + 62) + q(^i, £2 - 62) - 2q(sci, x 2 )) + . . . (7£ 
zo 2 



with 



+00 



G = Y G(x 2 + m 2 b 2 ) = -(F- 1 Q)(Q) 

7712 = — OO 
+ OO 

Gi= Y rn 2 b 2 G(x 2 + m 2 b 2 ) =i(F- 1 Q)'(0) 

7772 = — OO 
+ OO 

G 2 = E (m 2 b 2 ) 2 G(x 2 + m 2 b 2 ) = (F- l Q)"(0) (79) 

7772= — OO 

Relation (178!) involves a finite number of nodes around the point where the 
relation is written. It depends on the number of nodes chosen to approximate 
the boundary condition. This number can be 1 for a crude approximation in- 
volving only one node or can be larger. For a very large number of nodes, the 
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condition tends towards the true absorbing condition for a half-plane in the 
periodic media given by (176|) . Up to now everything has been written for peri- 
odic media but it is clear that homogeneous media are also periodic media and 
so all that has been said applies also to homogeneous media. The condition 
(1781) can be seen as a generalization of the Taylor approximation boundary 
condition proposed by [26J. But, while the boundary conditions in [26J were 
obtained by approximation of the exact continuous relations for specific prob- 
lems, they are obtained here directly and with general applicability from the 
discretized equations. 



3 Simple examples 



3. 1 Estimation of the accuracy 



In this section we try to estimate the quality of the proposed boundary condi- 
tion compared with known relations for the simple case of the two-dimensional 
acoustics. Consider first a plane wave incident on the plane y = at an angle 
9 with the normal to the plane. Let us define points at a horizontal distance D 
from the origin and with a vertical spacing h, see figure The sound pressure 
at point (D, Ih) is given by 

pi ^iK (cos 9D+sm9lh) (80) 



where K — u/c is the wavenumber and c is the sound velocity. The analytical 
force at the same points is the normal derivative in direction 1 given by 

f a = iK cos 9e iK{coseD+sineih) (81) 

For a point source at origin, the pressure is solution of 

Ap + K 2 p = -5(r) (82) 



where r is the distance from the origin and 5(.) is the Dirac function. The 
solution of this equation for the time dependence e~ lult is given by 

G(r) = l -H (Kr) (83) 
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where Hq is the Hankel function of zero order and first type. The analytical 
solution at each point I is 



Pa = iMK^D* + (lh) 2 ) 



34) 



and the analytical force at the same points is the normal derivative (in direc- 
tion 1) given by 



f 

J a 



iKD 



4JD 2 + (lh) 2 



H^KJD 2 + {lh)' 



55) 



The absorbing boundary condition described in the precedent section will 
allow to compute numerical forces f n at a node from the knowledge of p a . If 
the boundary condition was perfect one would have f l n = f l a but the proposed 
condition is approximate and one only has fn~f l a - The error can be estimated 
by 



I f 1 - f 1 1 

I J n J a I 



(86) 



The next step is to compute f l n from the method proposed in this paper and 
the error by relation ( J86l) to estimate the quality of the absorbing condition. 



3.2 Acoustic element 



Consider the rectangular four nodes acoustic element of size b\ x 6 2 - The 
elementary stiffness and mass matrices are given by 



K 



66i6 2 



(26| + 26?) (-26^ + 6?) 
(-261 + bj) (2bj + 26?) 
(-bt-bl) (6|-26?) 
(6 2 2 - 26?) (-bj-bl) 



[-bi 



(62 - 26?) 
(2b 2 + 26?) 



-26| 



(61 - 26?) 



-bl 



6?) (2b 2 



bl) 
2b\ + 6?) 
26?) 



37) 



M 



6i6 2 
36c 2 



4 2 12 
2 4 2 1 
12 4 2 
2 12 4 



(88) 
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and the dynamic stiffness matrix can then be determined by D = K — u; 2 M. 



It can be noted that the reduced set of displacement dofs q r contains only 
(\lb- Then, the matrices W and Wi have the following forms: 



W r 



1 e 



Wi 



1 e 







*9) 



The terms Aj in equation ([55]) are given by 



A (k)=Wl(k)B l (k)W (k) 
1 



12^ - 6^ + 2K%b 2 + ( 6^ + 6^ + KXh) cos(kb 2 ) 



A 1 (k)=W* (k)B l (k)W (k) 
1 



9 



- 6 T- + 2^ 2 &i6 2 + ( -3^ + 6^ + K 2 ^ ) cos(A;6 2 ) 
h b 2 V h b 2 I ~ 



A 2 (k)=W* 1 (k)B l (k)W 1 (k) 
1 



9 

= A,(k) 
A 3 (k) = W*(k)B l (k)W 1 {k) 
1 



- 6 T- + 2^ 2 6i6 2 + ( -3^ + 6^ + K 2 6!6 2 ) cos(A;6 2 ) 
b\ b 2 V oi »2 / 



1£ 

:A (A;) 



12^ - 6^ + 2^ 2 6 x 6 2 + (6^ + 6^ + K^bA cos(kb 2 ) 
b\ b 2 V b\ b 2 I 



(90) 



The eigensolutions of the spectral problem ([52]) are then determined and are 
given by 



in 



-i/i 



1 



2A 3 
1 

2A, 



-(A 1 +A 2 )±J(A 1 + A 2 y-4A 2 



-{A 1 + A 2 )T^{A l + A 2 f-AAl 



(91) 
(92) 



The signs are selected as in section 2.3. As there is only one dof in this case, 
one has n + = 1 and after normalization one can choose qi(fc) = 1. From 
relation (160]) . taking also Pj(— k) = 1, yields 



y(-k) 



A 2 + e'Mg 



As 
2A 3 



(-(A, + A 2 ) ± y j{A 1 + A 2 f - AAf) + A 2 1 
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±\J(A l + A 2 y-AAl 1 



this gives, with the notation of relation flTUj) . 



Q(-k) = ±-j(A 1 + A 2 )*-4A 2 

F(-k) = l (93) 



Near k = (this means in fact near the normal incidence) one has the devel- 
opment 

Mk) = -r - \^b x b 2 + 1 + 6^ + K%b 2 ) {kb 2 f + 0((kb 2 ) 4 ) 
b\ o Sb \ bi b 2 I 

Al{k) = \~ \ K% ^ + Yg (~ 3 ^ + 6 | + R2blb2 ) {kb2)2 + °^ kb ^ 

A 2 (k) = A 1 {k) 
A 3 (k) = A (k) 

and this leads to 



Q(-k) = ±zKb 2] jl-^(Kb i yx 

i - - 2 3 Yzmi — —h > + °^ kb2)A) (95) 



From relation (J58l) one also has 



r(Ar) = WSD, (Wo + e^Wj) 
= A 1 + e i,t A a 

= Q{-k) (96) 

Following the rule that the positive waves are such that Re {icuqj^fj} > 0, 
one has to choose the minus sign in relation (1931) . For the case Kbi <C 1 and 
k ^ K, one has the approximation: 

/(0) ~ -iKb 2 (97) 
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The power across the boundary is thus 

P = ^Re(f left ^ right v*) = ^Re(f r (0)(-iu)*) = l -Kb 2 u > (98) 

For the second order approximation one has 



G = -(F- 1 Q)(0) 
G 1 =i(F- 1 Q)'(0) 

G 2 = (F~ 1 Q)"(0) 

The relation between forces and displacements dofs on the boundary of the 
element is thus given by using (1751) 

f(x 1 ,x 2 )^iKb 2 q(x 1 ,x 2 ) 

+771^(9(^1) x 2 + h) + q(x 1 ,x 2 - b 2 ) - 2q(x u x 2 )) (100) 
2Kb 2 

At order one finds the classical approximation of the radiating boundary 
condition. The factor b 2 is present because the force is calculated over an edge 
of an element of length b 2 . 

We compare four solutions in the following 

(1) The zero order solution with the numerical computation of Q(0) leading 
to the relation f n = Gop l a 

(2) The zero order solution with the simplified computation of Q(0) leading 
to the relation f l n = iKb 2 p l a 

(3) The second order solution with the numerical computation of G 2 = 
(F -1 Q)"(0) w ((F- X Q)(5) + (F- x Q)(-5) - 2(F- 1 Q)(0))/5 2 leading to 

the relation f n = G oP l a + ^(p l a +1 + fc 1 - 2p l a ) 

(4) The second order solution with the simplify computation of Q(0) and 
(F~ 1 Q)"(0) leading to the relation f n = iKb 2 p l a + ^(p l + l +p l ~ 1 - 2p l a ) 

3. 3 Example 

Here we compute the error of relation (155]) for different cases as shown in figure 
El Case a) is for a sound pressure created by a plane wave while in case b) the 



iKb, 



12 



(Kb, 



iKbo 



ib 2 1 H 3 e— 



(AT 2 6i62) 2 
36 



K 



1 - ±{KK 



ib 2 
~K 



(99) 
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sound pressure is created by a point source. The acoustic element used for the 
computation of the boundary condition can be of size b\ x 6 2 = 0.01m x 0.01m 
or b\xb 2 = 0.05m x 0.05m. The sound velocity is c = 340m/s and the distance 
between the points is h = b 2 - 

The first example is for a pressure created by a plane wave at the incidence 
angle 9 = 10°. The error for the four cases listed in the precedent section are 
plotted in figure [6l The error is the same for any point on the vertical axis. It 
can be observed that the second order relations are much better than the first 
ones as expected. The comparison of the two sizes for the acoustic element 
shows that the size 0.05m x 0.05m can reduce the accuracy of the solution for 
high frequencies and the second order condition. In these cases it is better to 
use elements with small sizes. 

In figure [7] the error is plotted versus the angle of incidence of the plane wave. 
An acoustic element of size b\ x 6 2 = 0.01m x 0.01m has been used. The 
solution is accurate (error less than 1%) for angles up to 10° for a zero order 
condition and up to 30° for a second order condition. This clearly shows that 
second order conditions are much better for waves at oblique incidence. 

Finally the error is plotted versus the distance along the y axis in figure [8] for 
a pressure created by a point source at distance D = lm on the x axis and 
at the frequency 1000Hz. The reduction in accuracy can be observed as we 
move along the y axis leading to greater incidence angles in agreement with 
the precedent observation on the plane waves. All these points confirm the 
accuracy of the method proposed here. 



4 Finite element examples 

4-1 Acoustics 

In this section we use the precedent boundary condition to solve some finite 
element problems for different frequencies and mesh densities. We consider 
first a finite element acoustic problem on a square domain with a point source 
excitation in its center. A two dimensional domain of size lm x lm is generated 
by Ansys. The domain and an example of mesh are presented in figure M The 
size of the acoustic element is b = 0.025m leading to 40 x 40 elements for 
the whole domain. In each case, only square elements are used for the mesh. 
The sound velocity is c = 340m/s. Only mesh information, stiffness and mass 
matrices are picked out and are then introduced into Matlab to get the results 
presented below. The procedure is done over the frequency band [0, 2000Hz]. 
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Numerical Green's functions are calculated for zero and second order boundary 
conditions. The excitation is at point (0, 0) and the analytical solution in 
infinite space is given by formula (1841) . The Green's functions are presented 
in figure [TU] for a point at (0.3, 0) on the horizontal axis and in figure [TT] 
for a point at (0.3, 0.3) along the diagonal. Good agreements between the 
two types of absorbing boundary conditions and the analytical solution can 
be observed. Both boundary conditions fail at low frequencies because the 
size of the domain is too small compared to the wavelength. Similarly the 
error for high frequencies are of same orders for both boundary conditions. 
For intermediate frequencies the error is lower for the second order boundary 
condition. This is more clearly seen in figure [12] where the relative error for 
the point (0.3, 0.3) is plotted versus the frequency. 

The same results are presented in figure [J3] for the point (0.3, 0.3) and a 
mesh density of 80 x 80 elements. The solution is clearly much better at high 
frequencies meaning that the errors seen in figures [10] and [11] can be explained 
by elements too large for these frequencies and not by the quality of the 
boundary conditions. In figure [TU the domain is now 2m x 2m with 80 x 80 
elements, so with elements of the same size as for figures [10] and [11] Results 
are plotted for the point (0.3,0.3). Now the improvement is clearly seen for 
low frequencies while there is no difference for high frequencies. 



4-2 Two dimensional elastodynamics 



The analytical solution in direction e n of the two dimensional elastodynamics 
case when submitted to a unit force at origin in direction e m , is given by 



Gr. 



I 

4// 



Ad~ nm + B 



•En-Em 



(101) 



with: 



A = H (K T r) 
B-- 



1 



K T r 

2A+\H (K T r)+f3 2 H (K L r) 



where H and Hi are the Hankel functions of first type, of orders zero and 
one respectively. The wavenumbers are K L = uj/cl and K T = uj/ct for the 
longitudinal and transverse waves respectively. The velocities cl, ct and the 
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ratio between them (3 are given by, 

9 A + 2u , , Ev , -E 
cf = with A = - and u = — (102) 

_2 ^ 



In this example, the same global meshes as for the acoustic case are used. The 
boundary condition is computed from the square four nodes elements. The 
material is steel with E = 2.10 11 Pa, v = 0.3 and p = 7800/cg/m 3 and plane 
strain conditions are used in the computation. The sizes of the periodic cell 
can be b = 0.025m, 0.0125m or 0.00625m. 

In figure [151 numerical solutions are compared with analytical solutions for the 
point (0.5, 0) and (0, 0.5) for different sizes of the cell. The curves represent 
the real and imaginary parts of the first component of the displacement for an 
excitation at origin in direction 1. The same remarks as the previous examples 
can be made: 1) a denser mesh leads to lower errors over the high frequency 
band [8.10 3 Hz, 20.l0 3 Hz}; 2) for the low frequency band [0, 8.10 3 Hz], numer- 
ical results are different from the analytical solutions due to the finite size of 
the domain (L = lm). 

In figure [T6], the results for these two points are presented when the size of 
the domain is increased successively with L = lm, 2m and 4m. In this case, 
when the size of the cell is fixed to b = 0.025m, a larger domain leads to lower 
errors over the low frequency band [0, 8.10 3 Hz}. The results are the same as 
previously in the high frequency band because in this domain the precision 
depends on the size of the elements and not on the size of the global domain. 
Some improvements are however seen for intermediate frequencies. 

In figure [17] the error for boundary conditions of zero and second orders are 
plotted versus the frequency. The second order condition is much more accu- 
rate for intermediate frequencies as for the acoustic case. 



5 Conclusion 



In this paper, a method to determine absorbing boundary conditions for 
two dimensional periodic media has been presented. It works directly on the 
discretized equations. The boundary condition is first obtained as a global 
impedance relation and is then localized into boundary conditions of various 
orders. In two examples, good agreements were observed when compared with 
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analytical solutions. 



In any case, the proposed method is efficient because it requires only the dis- 
crete dynamic matrices which can be obtained by any standard finite element 
software. This method could be used for media with more complex behaviors 
than those presented in the precedent examples. 
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Fig. 2. Periodic medium near the exterior boundary. 
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Fig. 4. A cell in the periodic medium 



34 




Source 

o — 



case b) 



D 



x 
x 

X 
X 
X 
X 
X 



X 
X 
X 
X 



. 5. Points used to estimate the boundary condition with sound pressures created 
a plane wave (case a) and a point source (case b) 
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Fig. 6. Error versus the frequency for a plane wave at 10° for (a) an element size 
0.01m x 0.01m and (b) an element size 0.05m x 0.05m 
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Fig. 8. Error versus the distance for a point source 
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Configuration of the 2D domain 
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9. Example of finite element domain with a point source in its center. 
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Fig. 10. Comparison of analytical and numerical Green's functions, with the 

present method at order and at order 2 . . for the 2D acoustics at point 

(0.3, 0) with L = lm and b = 0.025m: a) real part, b) imaginary part. 
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a) 




Fig. 11. Comparison of analytical and numerical Green's functions, with the 

present method at order and at order 2 . . for the 2D acoustics at point 

(0.3, 0.3) with L = lm and b = 0.025m: a) real part, b) imaginary part. 
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Fig. 12. Comparison of relative errors for order and second order boundary 

conditions for the 2D acoustics at point (0.3, 0.3) with L = lm and b = 0.025m. 
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a) 




Fig. 13. Comparison of analytical and numerical Green's functions, with the 

present method at order and at order 2 . . for the 2D acoustics at point 

(0.3, 0.3) with L = lm and b = 0.0125m: a) real part, b) imaginary part. 
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Fig. 14. Comparison of analytical and numerical Green's functions, with the 

present method at order and at order 2 . . for the 2D acoustics at point 

(0.3, 0.3) with L = 2m and b = 0.025m: a) real part, b) imaginary part. 
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Fig. 15. Comparison of analytical and numerical Green's functions for the 2D elas- 
ticity with different sizes of elements: a) Real part at (0,0.5), b) Imaginary part at 
(0,0.5), c) Real part at (0.5,0), d) Imaginary part at (0.5,0). 
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Fig. 16. Comparison of analytical and numerical Green's functions for the 2D elas- 
ticity with different sizes of the domain: a) Real part at (0,0.5), b) Imaginary part 
at (0,0.5), c) Real part at (0.5,0), d) Imaginary part at (0.5,0). 



46 




Fig. 17. Comparison of relative errors for order and second order boundary 

conditions for the 2D elastodynamics at point (0.5, 0) with L = lm and b = 0.025m. 
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